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An  Algorithm  for  F-(y)  Using  Cubic  B-Splines 

O 

0.  Introduction 

This  note  describes  an  algorithm  whereby  the  distribution  of  the 
Maximum  of  the  Stationary  Gaussian  Markov  Process  over  an  interval  may  be 
computed  efficiently.  It  is  an  extension  of  the  earlier  report  (Keilson 
and  Ross,  1978  [1])  whose  notation  it  employs.  The  (zeros  and  residues) 
algorithm  of  the  earlier  report  is  one  of  the  starting  points  for  the 
development  of  the  new  algorithm.  The  relationship  of  the  old  and  new 
algorithms  is  described  in  the  first  two  sections  of  this  note.  Section  3 
provides  some  of  the  methodology  of  cubic  splines  and  Sections  4  and  5 
describe  its  use  in  the  algorithm.  Section  6  contains  some  comments  on 
the  accuracy  and  economy  of  the  method. 
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1 .  Problem  Statement 

The  function  Fa(y)  can  be  computed  to  better  than  6  decimal  place 

V 

accuracy  over  the  half  plane  0  <  0  <  ®,  -»  <  y  <  »  using  a  combination 
of  algorithms.  The  principal  ones  are  the  zeros  and  residues  method  for 
6  >  1  and  series  expansions  in  for  0  <  1.  Both  these  methods  require 
such  lengthy  computations  for  each  value  of  F„(y)  that  they  do  not  provide 
a  practical  algorithm  in  cases  where  considerable  numbers  of  values  of 
F  (y)  are  needed  cheaply  and  quicklv.  Nonetheless.  Fn(v)  is  a  simple 
and  well-behaved  function  of  its  areuments.  It  is  monotonicallv  increas¬ 
ing  as  a  function  of  y  for  all  0  and  monotonical ly  decreasing  as  a  function 
of  0  for  all  y.  It  therefore  seems  appropriate  to  look  for  some  represen¬ 
tation  of  the  function  which  is  easily  evaluated  and  does  not  require  too 
many  stored  constants.  The  latter  requirement  rules  out  the  straightforward 
idea  of  simply  storing  a  large  table  of  values  and  interpolating.  Use  may 
be  made  of  the  fact  that  over  much  of  the  0-y  plane  the  function  is  ade¬ 
quately  represented  by  a  simple  exponential  form,  but  since  in  some  regions 
this  does  not  give  even  one  decimal  place  accuracy,  additions  to  the  single 
exponential  form  are  required.  Since  the  principal  goals  for  this  algorithm 
are  low  cost,  speed  and  portability  rather  than  high  accuracy,  it  was  designed 
to  be  accurate  to  only  4  decimal  places  rather  than  the  6  places  achieved  by 
the  earlier,  slower  methods. 
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2.  Algorithm  Development 

Numerical  calculation  using  the  slow  version  of  the  algorithm  shows 
the  following  facts: 

(a)  For  y  <  -4,  F0(y)  <  0.0001  . 

(b)  For  y  >  -4  and  outside  the  rectangle  0  <  9  <  4  and  -4  <  y  <  4, 

F  (y)  may  be  represented  to  4  decimal  place  accuracy  by  the 
expression  ot(y)exp(-A (y)9) ,  where  A(y)  is  the  first  zero  in  the 
zeros  and  residues  representation  of  F0 (y)  and  a(v)  =  exp(-0.S  y*)* 
B(y)/(/2iT  A(y)).  6 Cy)  is  the  residue  of  D  j(-y)/D  ( -y)  at 

s  =  -  A (y)  (c.f .  equations  (2.11)  and  (2.10)  of  the  blue  1978 
report) . 

(c)  For  y  >  8,  a(y)  =  1.0000  and  A(y)  is  given  to  better  than  4  decimal 
place  accuracy  by 

y  2/2 

A (y)  =f  (1//27)  /  dx  ex 
0 

2 

=  (l//?n)e"[y  /2V/(1  +  1/y2  +  3/y4  +  15/y6  +  105/y8)  . 

(d)  Within  the  rectangle  0  <  9  <  4  and  -4  <  y  <  4,  F  (y)  may  be  written 

0 

as 


F0(y)  =  a(y)exp(-A(y)6)  +  G0(y) 

where  the  first  term  is  as  in  (b)  above,  and  G. (y)  is  a  correction 

0 

term  with  the  following  properties: 

I.  G0(y)  >  0 

II.  G0(y)  has  a  maximum  value  of  less  than  0.2  at  a  point 
close  to  9  =  0,  y  =  0.  It  falls  away  from  that  maximum 
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in  all  directions  and  in  less  than  0.0001  on  and 
outside  the  three  lines  y=-4,  0  <  H  4;  y  =  4, 

0  ^  9  <  4;  and  6=4,  -4  :£  y  £  4.  T 

In  view  of  facts  (a),  (b) ,  (c)  and  (d)  above,  the  problem  of  finding  a 
concise  representation  of  F Q(y)  becomes  one  of  finding: 

I,  a  representation  for  X(y)  and  ct(y)  for  -4  <  y  <  8 

II,  a  representation  for  G0(y)  throughout  the  area  -4<y<4,  O<0<4. 

In  both  cases  representation  in  terms  of  cubic  B-splines  are  used.  The  next 
section  is  a  description  of  the  properties  of  cubic  B-splines  needed  in  the 
development  of  the  algorithm. 


The  material  on  B-splines  used  in  the  development  of  the  algorithms 

may  be  found  in  Cox  [2],  Hayes  [3]  and  Hayes  and  Halliday  [4].  This  section 

summarizes  some  of  the  material  from  [2],  [3]  and  [4], 

Cubic  splines  provide  a  way  of  interpolating  a  function  given  at  a 

set  of  data  points  using  a  set  of  cubic  polynomial  arcs.  More  formally, 

a  cubic  spline  s(x)  on  a  set  of  knots  (t.  <  t.  <  ...  <  x  )  is 

r  1  2  n  1  2  n 

a  function  of  x  possessing  the  following  two  properties. 

I.  In  each  of  the  intervals  x  s  t, ;  t.  ,  <  x  <  x.,  j  =  2,...,n; 

1  3-1  V 

t^  <  x,  s(x)  is  a  polynomial  of  degree  3  or  lower. 

II.  S(x)  and  its  first  and  second  derivatives  are  continuous.  For  a 
given  finite  set  of  knots  the  set  of  cubic  splines  is  a  finite  dimensional 
linear  space.  In  the  development  of  algorithms  it  is  useful  to  work  with  a 
basis  set  of  splines  which  have  a  specific  form  and  which  lead  to  an  effi¬ 
cient  and  well-conditioned  evaluation  algorithm.  Then  anv  cubic  spline  on 

k  ’ 

that  knot  set  has  a  representation  of  the  form  S(x)  =  \  y.B.fxl  where 

k  i  =  l 

the  set  {B.}.  .  is  the  basis  set  and  y.  are  weights  chosen  to  make  S(x) 
ii=l  i  ° 

approximate  some  desired  function  f(x). 

A  commonly  used  basis  set  is  the  set  of  B-splines  which  may  be  defined 
as  follows  (cf.  [2],  [3],  [4]  and  [5]).  The  cubic  B-spline  B^(x)  is  the 
cubic  spline  with  knots  x.  ,,  x.  x.  x.  ,  and  x.,  which  is  zero  every- 

r  1-4  i-3  i-2  l-l  l  1 

where  except  in  the  range  x^  ^  <  x  <  x^.  This  defines  B^(x)  uniquely  except 
for  an  arbitrary  scale  factor  which  is  conventionally  chosen  to  make 

ao 

J  B^(x)dx  =  j.  A  cubic  B-spline  looks  like  a  piecewise  cubic  localized 
-00 

'hump'  function  as  shown  in  the  din'-ram  below. 
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The  function  and  its  first  two  derivatives  are  zero  at  the  end  points. 

Though  expressions  for  B^fx)  may  be  written  out  explicitly  in  terms  of 

powers  of  x  (cf.  [4],  p.  95),  these  do  not  always  provide  the  best  method 

for  evaluating  B^(x)  since  the  calculation  may  be  ill-conditioned  if  the 

knot  spacing  is  irregular.  A  stable  method  of  evaluating  B^(x)  uses 

recurrence  on  the  order  of  the  splines.  Splines  of  order  n+1  are  made  up 

of  polynomial  arcs  of  order  n,  so  cubic  splines  are  splines  of  order  4. 

Thus,  what  has  been  written  as  B^(x)  should  be  written  more  fully  as 

B  . (x) .  In  this  notation  the  recurrence  formula  to  be  used  (cf.  [2], 

^  >  ■* 

p.  137,  or  [3],  p.  149,  or  [4],  p.  95)  is 


B  .(x)  = 
n,  x 


(x  -  t  . )  B  .  .  (x)  +  (t.  -  x)  B  ,  .  (x) 
_ _ n-i  n-l,i-lv  i  n-l,i 

T  .  -  T  . 

i  i-n 


1/{Ti  -  Ti-P  if  Ti-1  s  x  "  Ti 


with 


0 


otherwise  . 
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Because  of  the  local  nature  of  B-splines  (see  diagram  above),  for  any 

given  x  no  more  than  4  of  the  members  of  the  basis  set  are  non- zero.  If 

x  is  a  knot,  only  3  are  non-zero.  Therefore,  to  evaluate  the  cubic  spline 
k 

S(x)  =  £  y.B.(x)  at  x  =  t  evaluate  the  4  B-splines  B. (t)  =  B.  .  (t)  ,  for 

i=l  1  1  1 

i  =  j , j+1 , j*2 , j  +  3  which  are  non-zero,  using  the  recurrence  at  the  top  of 

j  +  3 

the  page.  Then  s(t)  =  £  y.B.(t).  In  practice  it  is  convenient  to  eval- 

i=j  1  1 

uate  the  4  required  B-splines  together  via  a  calculation  scheme  which  looks 
like 


Step  1 

Step  2 

Step  3 

Step  4 

*i.j 

8-’n 

B3,  j 

B4 ,  j 

B2.M 

B3,j*l 

B4,j*l 

B7  .  . 
3,1*2 

B4 , j  +  2 

B4,j+3 

In  the  above  description  of  how  to  evaluate  a  cubic  spline,  it  has  beer. 

assumed  that  the  set  of  weights  {y.}^  has  been  given.  In  using  the  splii.e 
k  1  1  =  1 

S(x)  =  £  y.B.(x)  to  approximate  a  function  f(x),  an  appropriate  set  of 

i  =  l  1  1 

weights  must  be  calculated.  There  are  a  number  of  ways  to  do  that.  The  one 

chosen  depends  on  what  is  to  constitute  a  good  approximation  of  S(x)  to  f(\). 

For  the  present  purpose,  the  following  method  is  used.  Suppose  the  function 

f(x)  is  to  be  approximated  by  the  spline  s(x)  in  the  interval  a  <  x  <  b. 

Choose  n*6  knots  so  that  t.  <  t,  <  t.  <  t  ,  =  a  <  t.  ...  <  i  ,  = 

l  i=l  1  2  o  4  5  n+3 

b<t  .st  r  <  t  ,.  That  is,  n  of  the  knots  are  within  the  set  a  s  x  s  h 
n+4  n+5  n+6 

and  6  are  outside.  No  values  of  f(x)  are  needed  for  these  6  knots,  however. 

n+2 

Then  the  desired  spline  approximation  to  f(x)  is  s(x)  =  ][  y.B.(x)  where 

i  =  l  11 
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the  are  determined  by  the  n+2  condition  f(t.)  *  s(t^)  for  i  =  4,5, . . . ,n+3. 

It  t  I 

f  (t^)  =  s  (t4)  and  f  (t^^)  =  s  which  expresses  the  fact  that  the 

spline  matches  the  function  at  the  n  knots  in  a  £  x  s  b,  and  the  derivative 
of  the  spline  matches  the  derivative  of  the  function  at  the  end  points. 

Since  at  any  knot  only  3  of  the  B-splines  B^(x)  are  non-zero,  the  following 
set  of  equations  results: 


»  »  II 


YiW 

+  Y2B2(t4) 

+  Y3B3(T4) 

=  f  (t4) 

YiW 

+  WV 

+  WV 

= 

Y2B2(TS3 

+  Y3B3(T5} 

+  WV 

n+3^  +  Yn+lBn+l ^Tn+3^ 

+  Yn+2Bn+2(tn+3) 

'  f<V3> 

n+3^  +  Yn+lBn+l ^Tn+3^ 

Yn+2Bn+2^Tn+3^ 

=  f'tV3’ 

which  may  be  written  in  matrix  form  as 


By  =  f  . 

This  equation  may  be  solved  for  the  vector  y  by  inverting  the  matrix  B 
using  a  standard  Gaussian  elimination  algorithm.  An  alternative  method 
exploits  the  fact  that  §  is  a  band  matrix  with  no  more  than  three  non- zero 
elements  in  any  row  and  is  diagonally  dominant  so  simpler  and  more  efficient 
algorithms  may  be  employed  with  profit.  (For  details  see  Ahlberg,  Nilson 
and  Walsh  [6],  pp.  12-14.)  De  Boor  [5],  p.  206  solves  the  same  problem 
for  more  general  types  of  banded  diagonally  dominant  matrices.  Thus,  using 


any  of  a  number  of  more  or  less  standard  procedures  the  vector  of  weights 

n+2  n+2 

{ Y i } i_ ^  for  the  spline  approximation  S(x)  =  }  y.B.(x)  to  the  function 

1  i  =  l 

f(x)  may  be  computed. 

In  computing  a  spline  representation  the  optimal  choice  of  knot  posi¬ 
tions  which  gives  highest  accuracy  for  a  given  number  of  knots  is  quite 
complicated.  For  the  present  purpose,  equally  spaced  knots  were  used  and 
their  number  increased  (spacing  decreased)  until  the  desired  accuracy  was 
reached.  In  that  connection  it  may  be  noted  that  for  cubic  splines  and 
the  type  of  smooth  function  being  approximated,  the  maximum  error  in  the 
representation  is  proportional  to  the  fourth  power  of  the  know  separation. 
Thus,  increasing  the  number  of  knots  by  a  factor  of  2  improves  the  accuracy 
by  a  factor  of  16. 
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4.  Spline  representations  for  A(y),  a(y)  and  GQ(y) 

I.  A (y)  and  ct(y).  These  are  functions  of  the  single  variable  y  and 
the  application  of  the  algorithms  of  the  previous  section  is  straightforward. 
The  range  of  y  values  to  be  covered  was  -4  <  y  <  8.  Thirty-four  uniformly 
spaced  knots  were  used  (spacing  8/22  =  0.363636),  giving  36  spline  coeffi¬ 
cients.  Since  there  is  a  substantial  change  in  the  size  of  A (y)  and  a(y) 
over  the  range  of  y  (several  orders  of  magnitude  in  the  case  of  X(y)),  the 
functions  approximated  were  log(X(y))  and  log(a(y)).  The  values  of  the 
derivatives  at  the  end  points  were  approximated  using  a  5-point  numerical 
differentiation  formula. 

II.  G.(y)  is  a  function  of  2  variables.  For  each  of  23  equally  spaced 

values  of  y  (y  =  -4.0  +  8k/22;  k  =  0,1 ,...,22),  a  spline  representation  of 

G0(y)  as  a  function  of  /§"  in  the  range  0  <  /e  <  2  was  computed.  /§"  was 

used  rather  than  6  since  for  small  0,  Fg(y)  is  proportional  to  /e.  Eight 

knots  were  used  with  a  spacing  of  2/7  =  0.285714.  This  gave  10  spline 

coefficients,  each  of  which  is  a  function  of  y.  Each  of  these  functions  of 

y  were  then  spline  fitted  using  23  knots  giving  25  spline  coefficients.  On 

d  /J  2 

the  0=0  boundary  the  derivative  -  G  (y)  is  —  exp{-0.5y  },  which  is 

d/@" 

obtained  from  equation  (2.18)  on  page  14  of  the  blue  report  and  the  fact 

that  afy'1  exp(-A  (y)0)  =  0.  Derivatives  at  boundaries  other  than  0  =  0 
d/0 

were  computed  using  5  point  numerical  differentiation.  The  function  G„(y) 
is  thus  represented  by  an  array  of  250  coefficients  indexed  from 
i  =  1,  10  along  the  /e  dimension  and  from  j  =1,  25  along  the  y  dimension. 

The  algorithm  for  the  evaluation  of  G0(y)  follows  a  pattern  which  is 
the  reverse  of  what  has  just  been  described.  Four  spline  coefficients  are 


needed  to  calculate  G0(y)  as  a  function  of  /@,  and  they  are  obtained  by 

evaluating  their  spline  representation  as  a  function  of  y.  The  value  of 

G  (y)  is  therefore  obtained  from  a  4x4  block  of  coefficients  in  the  10x25 
b 

array  of  coefficients  {>-•}•  ,  ,  _  .  .  In  more  detail  the  procedure  is 

ij  i=l , 10; j=l , 25  * 

as  follows: 

I.  Find  K,  the  start  index  of  the  y  dimension  of  the  4x4  block  of 
spline  coefficients,  from  K  =  integer  part  of  2.75y  +  11.0. 

II.  Find  L,  the  start  index  of  the  0  dimension  of  the  4x4  block  of 
spline  coefficients,  from  L  =  integer  part  of  3.5/&. 

III.  Evaluate  the  4  spline  coefficients  in  the  spline  representation 
of  G0(y)  as  a  function  of  /@  for  a  given  value  of  v.  They  are 


K+3 


h  =  l  Y-.B  (y)  ,  j  =  L, 1+1 , L+2, L+3 
J  i  =  K  J 


L+3 


IV.  Evaluate  GQ(y)  =  l  <p  .B.(/e)  . 

j=L  J  1 


Some  economy  in  the  calculation  is  achieved  by  having  the  knot  set  in 
the  y  dimension  of  GQ(y)  coincide  with  the  knot  set  used  in  the  representa¬ 
tion  of  A(y)  and  a(y)  over  the  range  -4  <  y  <  4. 


-12- 


5 .  Later  modification  to  the  basic  algorithm 

After  the  algorithm  described  above  was  completed,  it  was  found  that 

the  accuracy  of  the  values  of  F0(y)  could  be  slightly  improved  by  enlarging 

the  domain  over  which  G.(y)  was  computed  by  one  knot  in  each  dimension.  The 

new  region  for  which  G  (y)  was  computed  was  thus  -4  <  y  s  4.363636, 

0 

0  <  /e  <  2.285714.  This  meant  that  the  number  of  spline  coefficients 
which  had  to  be  stored  increased  to  11x26  =  286.  At  the  same  time,  it 


was  realized  that  for  a  considerable  area  inside  that  region  GQ(y)  was  con¬ 
siderably  less  than  0.0001.  It  seemed  appropriate  to  simply  set  G  (v)  to 

zero  in  that  area  and  reduce  the  number  of  stored  constants  y. ..  The  number 

iJ 

of  constants  was  reduced  from  286  to  204,  and  the  actual  area  over  which 
G0(y)  is  computed  is  shown  in  diagram  1. 

The  penalty  paid  for  this  saving  of  storage  space  is  that  the  algorithm 
for  finding  the  block  of  16  constants  for  the  desired  y,  6  combination  is 
more  complicated.  The  storage  pattern  of  the  y^  no  longer  falls  into  a 
tidy  rectangular  grid.  Instead,  11  vectors  of  unequal  length  are  packed 
serially  in  a  one-dimensional  array.  Supplementary  arrays  are  needed  to 
mark  the  end  points  of  these  vectors.  It  is  questionable  whether  the  rela¬ 


tively  modest  saving  in  space  is  worth  the  added  complexity  incurred. 


DIAGRAM  OF  TME  HALF  PLANE  0  <  0  <  »,  -®  <  y  <  <*>  SHOWING  THE  REGIONS  INTO  WHICH  IT  IS  DIVIDED 
FOR  THF.  COMPUTATION  OF  Ffy)  BY  PROGRAM  FTHY .  REGIONS  EXTEND  TO  »  AT  TOP,  BOTTOM,  .AND  RIGHT 


6.  Accuracy  and  Economy 


The  knot  spacing  for  the  spline  representation  had  been  chosen  to 
make  the  algorithm  accurate  to  4  decimal  places.  Checks  were  run  which 
compared  the  values  given  by  the  spline-based  algorithm  with  the  earlier, 

6  decimal  place  accuracy  algorithm.  Comparisons  were  made  on  a  mesh  of 
points  at  y  =  -4.25(0.05)4.25,  0  =  0.0(0. 1)4. 5.  This  grid  is  considerably 
finer  than  the  knot  spacing.  Another  comparison  was  done  at  a  grid  made 
up  of  the  mid-points  between  knots  in  an  attempt  to  do  a  worst  case  com¬ 
parison.  In  both  these  comparisons  the  maximum  discrepancy  was  less  than 
0.000075. 

Since  cubic  splines  are  such  simple  functions  (cubics),  their  evalua¬ 
tion  is  rapid.  In  the  algorithm  8  B-splines  have  to  be  evaluated,  and 
these  can  be  done  in  2  blocks  of  4  using  the  recurrence  method.  The  bulk 
of  the  remainder  of  the  calculation  is  made  up  of  the  evaluation  of  2  expo¬ 
nentials,  one  square  root  and  7  inner  products,  each  with  4  components. 
Speed  tests  in  which  4000  representative  values  were  computed  showed  that 
the  algorithm  was  35  times  faster  than  the  earlier  6  decimal  place  version 
which  used  Chebychev  polynomial  approximation,  and  that  in  turn  was  about 
30  times  faster  than  the  original  zeros  and  residues  method.  Thus,  the 
spline  method  is  about  1000  times  faster  than  the  zeros  and  residues  method 
On  a  PDP11/34  with  hardware  floating  point  it  produces  about  200  values/sec 
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